Exploratory Analysis for SOAP by Sampling Plots
## by species types by PlotID
xtabs(~plotid+scientificname, data=SOAP.cat)
## scientificname
## plotid Abies concolor Amelanchier utahensis Arctostaphylos viscida
## SOAP1343 0 0 12
## SOAP139 0 0 14
## SOAP143 0 0 24
## SOAP1515 0 0 23
## SOAP1563 0 0 14
## SOAP1611 0 0 10
## SOAP1695 0 0 36
## SOAP187 0 0 0
## SOAP223 0 0 8
## SOAP255 0 0 15
## SOAP283 0 0 7
## SOAP299 0 0 3
## SOAP331 4 2 0
## SOAP43 7 0 0
## SOAP555 0 0 6
## SOAP63 0 0 6
## SOAP95 0 0 0
## SOAP991 0 0 1
## scientificname
## plotid Calocedrus decurrens Ceanothus cuneatus Ceanothus integerrimus
## SOAP1343 1 5 0
## SOAP139 72 0 1
## SOAP143 13 0 34
## SOAP1515 0 0 1
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 0 1 0
## SOAP187 81 0 0
## SOAP223 0 2 0
## SOAP255 14 0 0
## SOAP283 3 0 1
## SOAP299 3 0 0
## SOAP331 85 0 0
## SOAP43 63 0 0
## SOAP555 0 0 0
## SOAP63 38 0 0
## SOAP95 61 0 5
## SOAP991 1 41 0
## scientificname
## plotid Cercocarpus montanus var. glaber Corylus cornuta
## SOAP1343 2 0
## SOAP139 0 0
## SOAP143 0 0
## SOAP1515 0 0
## SOAP1563 0 0
## SOAP1611 0 0
## SOAP1695 0 0
## SOAP187 1 0
## SOAP223 16 0
## SOAP255 0 0
## SOAP283 0 0
## SOAP299 0 0
## SOAP331 0 0
## SOAP43 0 3
## SOAP555 0 0
## SOAP63 0 0
## SOAP95 0 0
## SOAP991 21 0
## scientificname
## plotid Pinus lambertiana Pinus ponderosa Prunus emarginata
## SOAP1343 0 1 0
## SOAP139 3 6 0
## SOAP143 0 38 0
## SOAP1515 0 1 0
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 0 4 0
## SOAP187 0 8 0
## SOAP223 0 0 0
## SOAP255 0 51 0
## SOAP283 0 0 0
## SOAP299 0 16 0
## SOAP331 1 11 0
## SOAP43 9 6 0
## SOAP555 0 0 0
## SOAP63 0 53 0
## SOAP95 3 18 0
## SOAP991 0 0 1
## scientificname
## plotid Quercus chrysolepis Quercus garryana Quercus kelloggii
## SOAP1343 0 1 1
## SOAP139 19 0 2
## SOAP143 0 0 1
## SOAP1515 1 0 0
## SOAP1563 0 0 0
## SOAP1611 0 0 0
## SOAP1695 50 0 0
## SOAP187 2 0 2
## SOAP223 14 1 1
## SOAP255 15 0 9
## SOAP283 3 0 0
## SOAP299 4 0 0
## SOAP331 6 0 23
## SOAP43 0 0 9
## SOAP555 0 0 0
## SOAP63 10 0 15
## SOAP95 17 0 2
## SOAP991 0 1 0
## scientificname
## plotid Rhamnus ilicifolia Ribes roezlii
## SOAP1343 0 1
## SOAP139 0 5
## SOAP143 0 4
## SOAP1515 0 0
## SOAP1563 0 0
## SOAP1611 0 0
## SOAP1695 12 3
## SOAP187 0 1
## SOAP223 0 1
## SOAP255 0 0
## SOAP283 0 0
## SOAP299 0 0
## SOAP331 0 0
## SOAP43 0 3
## SOAP555 0 0
## SOAP63 0 0
## SOAP95 0 0
## SOAP991 0 1
## aggregate summary stats
# Re-instate PlotID
SOAP.cont$plotid<-insitu.SOAP$plotid
# group by
melted <- melt(SOAP.cont, id.vars=c("plotid"))
SOAP.aggregate<-summarise(group_by(melted, plotid, variable),
mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE),
min=min(value, na.rm=TRUE), max=max(value, na.rm=TRUE))
SOAP.aggregate
## Source: local data frame [126 x 6]
## Groups: plotid [?]
##
## plotid variable mean sd min max
## (chr) (fctr) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 dbh 8.400000 6.437391 1.0 15.4
## 2 SOAP1343 dbhheight 120.000000 33.879582 10.0 130.0
## 3 SOAP1343 basalcanopydiam 29.287500 30.840844 0.0 140.0
## 4 SOAP1343 maxcanopydiam 3.920833 5.736684 0.6 30.0
## 5 SOAP1343 stemheight 2.112500 1.212368 0.7 6.5
## 6 SOAP1343 livingcanopy 71.625000 33.217089 0.0 100.0
## 7 SOAP1343 inplotcanopy 100.000000 0.000000 100.0 100.0
## 8 SOAP139 dbh 7.274510 19.207526 0.4 159.8
## 9 SOAP139 dbhheight 80.983607 59.089648 10.0 130.0
## 10 SOAP139 basalcanopydiam 5.311475 13.882674 0.0 60.0
## .. ... ... ... ... ... ...
## WHOA-- That's a big spicy meatball!
## Break it down by plot again-- this time stratify by parameter
#overside columnn width restriction
options(dplyr.width = Inf)
# dbh
insitu_dbh <- SOAP.cont %>%
group_by(plotid) %>%
summarise(dbh.min = min(dbh, na.rm=TRUE), dbh.max = max(dbh, na.rm=TRUE), dbh.mean=mean(dbh, na.rm=TRUE), dbh.median=median(dbh, na.rm=TRUE),dbh.sd=sd(dbh, na.rm=TRUE))
insitu_dbh
## Source: local data frame [18 x 6]
##
## plotid dbh.min dbh.max dbh.mean dbh.median dbh.sd
## (chr) (dbl) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 1.0 15.4 8.400000 8.60 6.437391
## 2 SOAP139 0.4 159.8 7.274510 1.70 19.207526
## 3 SOAP143 0.4 44.0 6.018889 1.40 11.503276
## 4 SOAP1515 1.0 24.0 11.646154 10.50 6.357753
## 5 SOAP1563 0.0 0.0 0.000000 0.00 0.000000
## 6 SOAP1611 0.0 0.0 0.000000 0.00 0.000000
## 7 SOAP1695 0.0 60.7 7.654667 1.90 12.417656
## 8 SOAP187 1.2 102.9 8.220652 2.75 16.507677
## 9 SOAP223 0.6 107.9 16.823529 10.20 26.028387
## 10 SOAP255 0.4 78.3 10.751724 7.30 13.063271
## 11 SOAP283 0.0 33.5 6.492857 0.00 12.273578
## 12 SOAP299 1.1 64.1 32.685000 33.90 21.016767
## 13 SOAP331 0.4 72.1 10.714729 5.90 12.610971
## 14 SOAP43 0.7 84.0 11.066316 5.80 16.032312
## 15 SOAP555 0.0 0.0 0.000000 0.00 0.000000
## 16 SOAP63 0.2 50.8 14.746364 13.25 12.035816
## 17 SOAP95 0.3 61.3 10.762626 5.50 12.969474
## 18 SOAP991 1.4 10.5 6.420000 8.20 4.490212
## functionalization might be best here...but for now
## JV scripting
# dbhheight
insitu_dbhheight <- SOAP.cont %>%
group_by(plotid) %>%
summarise(dbhheight.min = min(dbhheight, na.rm=TRUE), dbhheight.max = max(dbhheight, na.rm=TRUE), dbhheight.mean=mean(dbhheight, na.rm=TRUE), dbhheight.median=median(dbhheight, na.rm=TRUE),dbhheight.sd=sd(dbhheight, na.rm=TRUE))
insitu_dbhheight
## Source: local data frame [18 x 6]
##
## plotid dbhheight.min dbhheight.max dbhheight.mean dbhheight.median
## (chr) (dbl) (dbl) (dbl) (dbl)
## 1 SOAP1343 10.0 130 120.00000 130.0
## 2 SOAP139 10.0 130 80.98361 130.0
## 3 SOAP143 10.0 130 59.47368 10.0
## 4 SOAP1515 1.0 130 70.88462 130.0
## 5 SOAP1563 0.0 0 0.00000 0.0
## 6 SOAP1611 0.0 0 0.00000 0.0
## 7 SOAP1695 0.0 130 66.36604 38.5
## 8 SOAP187 2.0 130 69.07368 10.0
## 9 SOAP223 10.0 130 94.18605 130.0
## 10 SOAP255 10.0 140 109.32692 130.0
## 11 SOAP283 0.0 130 28.18571 0.0
## 12 SOAP299 10.0 130 125.38462 130.0
## 13 SOAP331 10.0 130 83.65909 130.0
## 14 SOAP43 10.0 130 103.60000 130.0
## 15 SOAP555 0.0 0 0.00000 0.0
## 16 SOAP63 0.8 130 111.23607 130.0
## 17 SOAP95 10.0 130 85.84906 130.0
## 18 SOAP991 10.0 130 121.04478 130.0
## dbhheight.sd
## (dbl)
## 1 33.87958
## 2 59.08965
## 3 59.33022
## 4 65.13944
## 5 0.00000
## 6 0.00000
## 7 61.59260
## 8 60.22498
## 9 55.12896
## 10 45.67221
## 11 55.19133
## 12 23.53394
## 13 58.62428
## 14 49.95998
## 15 0.00000
## 16 43.87607
## 17 58.14379
## 18 31.77260
# basalcanopydiam
insitu_basalcanopydiam <- SOAP.cont %>%
group_by(plotid) %>%
summarise(basalcanopydiam.min = min(basalcanopydiam, na.rm=TRUE), basalcanopydiam.max = max(basalcanopydiam, na.rm=TRUE), basalcanopydiam.mean=mean(basalcanopydiam, na.rm=TRUE), basalcanopydiam.median=median(basalcanopydiam, na.rm=TRUE),basalcanopydiam.sd=sd(basalcanopydiam, na.rm=TRUE))
insitu_basalcanopydiam
## Source: local data frame [18 x 6]
##
## plotid basalcanopydiam.min basalcanopydiam.max basalcanopydiam.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.0 140 29.2875000
## 2 SOAP139 0.0 60 5.3114754
## 3 SOAP143 0.0 40 2.3771930
## 4 SOAP1515 0.0 117 23.8461538
## 5 SOAP1563 2.0 90 27.2142857
## 6 SOAP1611 1.1 220 42.6100000
## 7 SOAP1695 0.0 94 10.3000000
## 8 SOAP187 0.0 48 1.0315789
## 9 SOAP223 0.0 58 11.3255814
## 10 SOAP255 0.0 80 5.5000000
## 11 SOAP283 0.0 31 14.0000000
## 12 SOAP299 0.0 125 15.9230769
## 13 SOAP331 0.0 60 0.5174242
## 14 SOAP43 0.0 38 0.8320000
## 15 SOAP555 7.0 130 90.8333333
## 16 SOAP63 0.0 61 2.3360656
## 17 SOAP95 0.0 30 0.7809524
## 18 SOAP991 0.0 163 34.6567164
## basalcanopydiam.median basalcanopydiam.sd
## (dbl) (dbl)
## 1 22.5 30.840844
## 2 0.0 13.882674
## 3 0.0 6.012362
## 4 20.0 21.502916
## 5 24.0 22.874839
## 6 18.0 64.730286
## 7 0.0 17.087255
## 8 0.0 6.104502
## 9 6.0 15.085030
## 10 0.0 15.270379
## 11 13.0 11.948608
## 12 0.0 35.305436
## 13 0.0 5.253389
## 14 0.0 4.755664
## 15 110.0 48.803347
## 16 0.0 8.759987
## 17 0.0 3.855518
## 18 24.0 31.471754
# maxcanopydiam
insitu_maxcanopydiam <- SOAP.cont %>%
group_by(plotid) %>%
summarise(maxcanopydiam.min = min(maxcanopydiam, na.rm=TRUE), maxcanopydiam.max = max(maxcanopydiam, na.rm=TRUE), maxcanopydiam.mean=mean(maxcanopydiam, na.rm=TRUE), maxcanopydiam.median=median(maxcanopydiam, na.rm=TRUE),maxcanopydiam.sd=sd(maxcanopydiam, na.rm=TRUE))
insitu_maxcanopydiam
## Source: local data frame [18 x 6]
##
## plotid maxcanopydiam.min maxcanopydiam.max maxcanopydiam.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.6 30.0 3.920833
## 2 SOAP139 0.2 30.0 1.752459
## 3 SOAP143 0.2 11.5 1.240351
## 4 SOAP1515 0.7 5.4 2.334615
## 5 SOAP1563 0.5 3.9 1.885714
## 6 SOAP1611 1.0 5.1 2.240000
## 7 SOAP1695 0.1 60.0 3.513208
## 8 SOAP187 0.0 80.0 2.309474
## 9 SOAP223 0.3 10.1 2.576744
## 10 SOAP255 0.0 12.6 3.031731
## 11 SOAP283 0.0 3.8 1.807143
## 12 SOAP299 0.0 10.1 4.330769
## 13 SOAP331 0.0 8.5 1.934351
## 14 SOAP43 0.0 18.1 2.440000
## 15 SOAP555 1.3 4.0 2.200000
## 16 SOAP63 0.0 6.8 2.604918
## 17 SOAP95 0.0 12.4 2.096226
## 18 SOAP991 0.4 4.7 1.962687
## maxcanopydiam.median maxcanopydiam.sd
## (dbl) (dbl)
## 1 2.85 5.7366843
## 2 0.90 3.2085684
## 3 0.60 1.8008709
## 4 2.10 1.0403624
## 5 1.80 1.0021954
## 6 1.80 1.3150412
## 7 2.20 7.1522634
## 8 0.80 8.3154186
## 9 1.90 2.1445209
## 10 2.40 2.3197449
## 11 1.70 1.2256642
## 12 4.25 2.9714332
## 13 1.40 1.6596928
## 14 1.95 2.4516743
## 15 1.95 0.9423375
## 16 2.45 1.6878688
## 17 0.90 2.2582304
## 18 2.00 0.9766727
# stemheight
insitu_stemheight <- SOAP.cont %>%
group_by(plotid) %>%
summarise(stemheight.min = min(stemheight, na.rm=TRUE), stemheight.max = max(stemheight, na.rm=TRUE), stemheight.mean=mean(stemheight, na.rm=TRUE), stemheight.median=median(stemheight, na.rm=TRUE),stemheight.sd=sd(stemheight, na.rm=TRUE))
insitu_stemheight
## Source: local data frame [18 x 6]
##
## plotid stemheight.min stemheight.max stemheight.mean
## (chr) (dbl) (dbl) (dbl)
## 1 SOAP1343 0.7 6.5 2.1125000
## 2 SOAP139 0.5 120.0 4.7459016
## 3 SOAP143 0.5 19.7 2.5307018
## 4 SOAP1515 0.6 3.4 1.7769231
## 5 SOAP1563 0.0 1.9 1.0714286
## 6 SOAP1611 0.6 4.2 1.7200000
## 7 SOAP1695 0.6 17.0 3.0566038
## 8 SOAP187 0.5 134.0 4.8789474
## 9 SOAP223 0.5 8.2 2.4441860
## 10 SOAP255 0.0 44.7 5.1019231
## 11 SOAP283 0.0 13.8 3.3285714
## 12 SOAP299 1.0 33.1 13.8961538
## 13 SOAP331 0.5 39.2 5.0098485
## 14 SOAP43 0.5 51.1 5.0470000
## 15 SOAP555 0.8 1.4 0.9666667
## 16 SOAP63 0.6 33.0 9.2327869
## 17 SOAP95 0.5 28.1 6.0509434
## 18 SOAP991 -1.6 3.8 1.8537313
## stemheight.median stemheight.sd
## (dbl) (dbl)
## 1 2.05 1.2123683
## 2 1.50 13.0803565
## 3 1.00 4.1863057
## 4 1.70 0.5867249
## 5 1.20 0.5426836
## 6 1.30 1.2318008
## 7 2.20 2.8116770
## 8 1.40 14.9537862
## 9 1.60 1.9816182
## 10 2.35 6.7463328
## 11 1.80 4.1795157
## 12 10.30 12.8098550
## 13 2.40 6.4031226
## 14 2.95 7.8090010
## 15 0.90 0.2250926
## 16 5.00 9.0170403
## 17 1.70 8.2020583
## 18 1.90 0.9241422
# livingcanopy
insitu_livingcanopy <- SOAP.cont %>%
group_by(plotid) %>%
summarise(livingcanopy.min = min(livingcanopy, na.rm=TRUE), livingcanopy.max = max(livingcanopy, na.rm=TRUE), livingcanopy.mean=mean(livingcanopy, na.rm=TRUE), livingcanopy.median=median(livingcanopy, na.rm=TRUE),livingcanopy.sd=sd(livingcanopy, na.rm=TRUE))
insitu_livingcanopy
## Source: local data frame [18 x 6]
##
## plotid livingcanopy.min livingcanopy.max livingcanopy.mean
## (chr) (int) (int) (dbl)
## 1 SOAP1343 0 100 71.62500
## 2 SOAP139 0 100 79.71311
## 3 SOAP143 0 100 95.92105
## 4 SOAP1515 20 100 82.42308
## 5 SOAP1563 NA NA NaN
## 6 SOAP1611 NA NA NaN
## 7 SOAP1695 0 100 88.57576
## 8 SOAP187 0 100 23.15789
## 9 SOAP223 0 100 79.48837
## 10 SOAP255 0 100 66.89320
## 11 SOAP283 NA NA NaN
## 12 SOAP299 0 100 84.61538
## 13 SOAP331 0 100 78.78788
## 14 SOAP43 0 100 81.90000
## 15 SOAP555 NA NA NaN
## 16 SOAP63 0 100 62.58197
## 17 SOAP95 0 100 84.90566
## 18 SOAP991 0 100 71.41791
## livingcanopy.median livingcanopy.sd
## (dbl) (dbl)
## 1 80 33.21709
## 2 100 39.42405
## 3 100 18.47703
## 4 90 20.29418
## 5 NA NaN
## 6 NA NaN
## 7 100 19.37172
## 8 0 42.40793
## 9 90 24.98225
## 10 100 46.67953
## 11 NA NaN
## 12 100 36.79465
## 13 100 41.03676
## 14 100 38.57814
## 15 NA NaN
## 16 100 48.22568
## 17 100 35.96944
## 18 80 26.73957
# inplotcanopy
insitu_inplotcanopy <- SOAP.cont %>%
group_by(plotid) %>%
summarise(inplotcanopy.min = min(inplotcanopy, na.rm=TRUE), inplotcanopy.max = max(inplotcanopy, na.rm=TRUE), inplotcanopy.mean=mean(inplotcanopy, na.rm=TRUE), inplotcanopy.median=median(inplotcanopy, na.rm=TRUE),inplotcanopy.sd=sd(inplotcanopy, na.rm=TRUE))
insitu_inplotcanopy
## Source: local data frame [18 x 6]
##
## plotid inplotcanopy.min inplotcanopy.max inplotcanopy.mean
## (chr) (int) (int) (dbl)
## 1 SOAP1343 100 100 100.00000
## 2 SOAP139 90 100 99.91803
## 3 SOAP143 100 100 100.00000
## 4 SOAP1515 100 100 100.00000
## 5 SOAP1563 NA NA NaN
## 6 SOAP1611 NA NA NaN
## 7 SOAP1695 65 100 98.03030
## 8 SOAP187 100 100 100.00000
## 9 SOAP223 100 100 100.00000
## 10 SOAP255 100 100 100.00000
## 11 SOAP283 NA NA NaN
## 12 SOAP299 100 100 100.00000
## 13 SOAP331 100 100 100.00000
## 14 SOAP43 40 100 98.90000
## 15 SOAP555 NA NA NaN
## 16 SOAP63 100 100 100.00000
## 17 SOAP95 100 100 100.00000
## 18 SOAP991 100 100 100.00000
## inplotcanopy.median inplotcanopy.sd
## (dbl) (dbl)
## 1 100 0.0000000
## 2 100 0.9053575
## 3 100 0.0000000
## 4 100 0.0000000
## 5 NA NaN
## 6 NA NaN
## 7 100 6.7867965
## 8 100 0.0000000
## 9 100 0.0000000
## 10 100 0.0000000
## 11 NA NaN
## 12 100 0.0000000
## 13 100 0.0000000
## 14 100 7.7713538
## 15 NA NaN
## 16 100 0.0000000
## 17 100 0.0000000
## 18 100 0.0000000
create a qualitative landscape assessment index (completely experimental- BE NICE!)
#subset SOAP.grand data set
#s<-subset(SOAP.grand, select = c(plot.id, easting, northing, max.dbh, max.dbhheight, max.basalcanopydiam, max.maxcanopydiam, max.stemheight))
## categorize by quartiles, use quantile values from summary and assign score
# dbh
summary(SOAP.grand$max.dbh)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 41.38 62.70 65.91 88.72 159.80
dbh.cat <- cut(SOAP.grand$max.dbh, breaks = c(-1, 41.38, 62.70, 88.72, 159.8))
dbh.cat<-unclass(dbh.cat)
#basalcanopydiam
summary(SOAP.grand$max.basalcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30.00 39.50 59.00 77.83 77.00 220.00
max.basalcanopydiam.cat<- cut(SOAP.grand$max.basalcanopydiam, breaks= c(0, 39.50, 59.0, 77.00, 222.00))
max.basalcanopydiam.cat<- unclass(max.basalcanopydiam.cat)
# maxcanopydiam
summary(SOAP.grand$max.maxcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 3.80 5.95 10.10 17.51 15.25 80.00 1
## impute missing value at Plot SOAP331
SOAP.grand$max.maxcanopydiam[8] <- mean(SOAP.grand$max.maxcanopydiam, na.rm=TRUE)
# determine class quantiles scores
summary(SOAP.grand$max.maxcanopydiam)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.800 6.375 10.800 17.510 17.660 80.000
max.maxcanopydiam.cat<-cut(SOAP.grand$max.maxcanopydiam, breaks= c(0, 6.375, 10.80, 17.660, 80.00))
max.maxcanopydiam.cat<-unclass(max.maxcanopydiam.cat)
#stem height
summary(SOAP.grand$max.stemheight)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.80 12.40 30.55 40.68 42.18 134.00
max.stemheight.cat<-cut(SOAP.grand$max.stemheight, breaks= c(0,12.40,30.55, 42.18, 134))
max.stemheight.cat<-unclass(max.stemheight.cat)
s<-data.frame(dbh.cat, max.basalcanopydiam.cat,max.basalcanopydiam.cat,max.stemheight.cat)
tree.score<-rowSums(s)
s$tree.score<-tree.score
s$tree.score.perc<- (tree.score/16)
## add on to SOAP.grand dataset
SOAP.grand<-c(SOAP.grand,s)
boxplot(s$tree.score.perc)
## simple linear regression using tree score ()
mod5.ndvi<-lm(tree.score.perc~SOAP_NDVI, data=SOAP.grand)
## also- makesure you append CHM values from other data stack to model
## import chm
SOAP.chm <- raster("../NEONdata/D17-California/SOAP/2013/lidar/SOAP_lidarCHM.tif")
SOAP.chm[SOAP.chm==0]<- NA
plot(SOAP.chm, main="Canopy Height Model for SOAP Site")
# Insitu sampling took place within 40m x 40m square plots, so we use a 20m radius.
# Note that below will return a dataframe containing the max height
# calculated from all pixels in the buffer for each plot (NEON Website)
## we will use these values to populate in-situ data set
chm.max <- extract(SOAP.chm,
SOAP_plots,
buffer = 20,
fun=(max),
sp=TRUE,
stringsAsFactors=FALSE)
# CHM distribution visualization for each plot
# created histograms for all 14 plots of data
cent_ovrList <- extract(SOAP.chm,SOAP_plots,buffer = 20)
for (i in 1:14) {
hist(cent_ovrList[[i]], main=(paste("plot",i)))
}
## merge data into data frame, can use approach to create spatial object
## But purposes of capstone will conduct "boring regression analysis"
## merge grouped data (insitu.max) and extracted NDVI (NDVI.max)
## using plotid to merge
## first must append NDVI.max file with SOAP
#create data frame for NDVI.max
chm.max.df<-chm.max@data
# add SOAP prefix to IDs
chm.max.df$plot.id<- "SOAP"
chm.max.df$plot.id<-(sapply(chm.max.df$plot.id, paste0, chm.max.df$ID))[,1]
# rename column
colnames(chm.max.df)[colnames(chm.max.df) == 'SOAP'] <- 'Plot.ID'
SOAP.grand<-merge(chm.max.df, SOAP.grand,
by.x = 'plot.id',
by.y = 'plot.id')
# F* Yes
## Response variables: dbh, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy
## Predictor variable: SOAP_lidarCHM
## Simple Linear Fit
mod.chm<-lm(max.dbh~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.dbh, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" DBH (cm)")
abline(mod.chm)
summary(mod.chm)
##
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.517 -16.469 -8.646 0.459 78.467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9555 30.4248 0.163 0.874
## SOAP_lidarCHM 1.8172 0.8449 2.151 0.057 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 38.35 on 10 degrees of freedom
## Multiple R-squared: 0.3163, Adjusted R-squared: 0.2479
## F-statistic: 4.626 on 1 and 10 DF, p-value: 0.05698
mod2.chm<-lm(max.basalcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.basalcanopydiam, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Base Canopy Diameter (cm)")
abline(mod2.chm)
summary(mod2.chm)
##
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -70.543 -47.516 -5.578 19.236 110.410
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 133.947 46.005 2.912 0.0155 *
## SOAP_lidarCHM -1.673 1.278 -1.309 0.2197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.99 on 10 degrees of freedom
## Multiple R-squared: 0.1464, Adjusted R-squared: 0.06102
## F-statistic: 1.715 on 1 and 10 DF, p-value: 0.2197
mod3.chm<-lm(max.maxcanopydiam~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.maxcanopydiam, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Maximum canopy diameter(m)")
abline(mod3.chm)
summary(mod3.chm)
##
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.240 -9.587 -3.554 0.542 49.986
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -11.5842 14.4213 -0.803 0.4405
## SOAP_lidarCHM 0.8674 0.4005 2.166 0.0556 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.18 on 10 degrees of freedom
## Multiple R-squared: 0.3193, Adjusted R-squared: 0.2512
## F-statistic: 4.691 on 1 and 10 DF, p-value: 0.05556
mod4.chm<-lm(max.stemheight~SOAP_lidarCHM, data=SOAP.grand)
plot(SOAP.grand$max.stemheight, SOAP.grand$SOAP_lidarCHM,
xlab="Max CHM", ylab=" Stem Height (m)")
abline(mod4.chm)
summary(mod4.chm)
##
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.832 -17.010 -5.808 9.709 59.479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -38.0415 23.7283 -1.603 0.13997
## SOAP_lidarCHM 2.3470 0.6589 3.562 0.00516 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.91 on 10 degrees of freedom
## Multiple R-squared: 0.5592, Adjusted R-squared: 0.5151
## F-statistic: 12.69 on 1 and 10 DF, p-value: 0.005165
## add on to SOAP.grand dataset
SOAP.grand<-c(SOAP.grand,s)
boxplot(s$tree.score.perc)
## simple linear regression using tree score
mod5.chm<-lm(tree.score.perc~SOAP_lidarCHM, data=SOAP.grand)
mod5.chm
##
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM, data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM
## 0.420212 0.006105
## Full model in respect to dbh
mod.full<- lm(max.dbh~SOAP_lidarCHM + SOAP_NDVI, data=SOAP.grand)
mod.full
##
## Call:
## lm(formula = max.dbh ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM SOAP_NDVI
## -137.057 1.505 172.771
## Full model in respect to base canopy diameter
mod2.full<-lm(max.basalcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod2.full
##
## Call:
## lm(formula = max.basalcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI,
## data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM SOAP_NDVI
## 856.50684 -0.08263 -879.05838
## Full model in respect to max canopy
mod3.full<-lm(max.maxcanopydiam~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod3.full
##
## Call:
## lm(formula = max.maxcanopydiam ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM SOAP_NDVI
## -15.1373 0.8595 4.3227
## Full model in respect to stem height
mod4.full<-lm(max.stemheight~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod4.full
##
## Call:
## lm(formula = max.stemheight ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM SOAP_NDVI
## -9.646 2.410 -34.546
## Full model in respect to % score
mod5.full<- lm(tree.score.perc~SOAP_lidarCHM +SOAP_NDVI, data=SOAP.grand)
mod5.full
##
## Call:
## lm(formula = tree.score.perc ~ SOAP_lidarCHM + SOAP_NDVI, data = SOAP.grand)
##
## Coefficients:
## (Intercept) SOAP_lidarCHM SOAP_NDVI
## 1.446148 0.008363 -1.248143
#descriptive statistics for in-situ field data products
stargazer(SOAP.cont, type = "text", title="SOAP-wide Descriptive statistics", digits=1, out="SOAPwideStats.txt")
##
## SOAP-wide Descriptive statistics
## ==============================================
## Statistic N Mean St. Dev. Min Max
## ----------------------------------------------
## dbh 982 10.1 14.9 0.0 159.8
## dbhheight 1,231 86.3 58.3 0.0 140.0
## basalcanopydiam 1,230 7.6 19.3 0.0 220.0
## maxcanopydiam 1,230 2.4 3.8 0.0 80.0
## stemheight 1,231 4.8 8.5 -1.6 134.0
## livingcanopy 1,146 74.2 41.1 0 100
## inplotcanopy 1,147 99.8 2.9 40 100
## ----------------------------------------------
stargazer(data.frame(insitu.max), type = "text", title="Maximum Value Plot-wide Descriptive statistics", digits=1, out="PlotwideStats-freq.txt")
##
## Maximum Value Plot-wide Descriptive statistics
## ===============================================
## Statistic N Mean St. Dev. Min Max
## -----------------------------------------------
## max.dbh 18 53.9 43.4 0.0 159.8
## max.dbhheight 18 108.9 50.2 0 140
## max.basalcanopydiam 18 88.1 52.0 30 220
## max.maxcanopydiam 17 18.1 21.4 3.8 80.0
## max.stemheight 18 31.3 38.3 1.4 134.0
## max.livingcanopy 14 100.0 0.0 100 100
## max.inplotcanopy 14 100.0 0.0 100 100
## -----------------------------------------------
## regression products
## NDVI-only
stargazer(mod.ndvi, mod2.ndvi, mod3.ndvi, mod4.ndvi, mod5.ndvi, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("NDVI"), ci=TRUE, note= "NDVI as the only predictor. Maximum NDVI values determined from 20 ft buffers at each respective SOAP station ", out="NDVImodels2.htm")
##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NDVI</td><td>394.748</td><td>-891.248<sup>***</sup></td><td>132.966</td><td>320.930</td><td>-0.014</td></tr>
## <tr><td style="text-align:left"></td><td>(-122.167, 911.663)</td><td>(-1,433.297, -349.200)</td><td>(-142.008, 407.940)</td><td>(-197.733, 839.592)</td><td>(-2.300, 2.272)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-282.516</td><td>864.495<sup>***</sup></td><td>-100.070</td><td>-242.586</td><td>0.638</td></tr>
## <tr><td style="text-align:left"></td><td>(-739.388, 174.355)</td><td>(385.409, 1,343.581)</td><td>(-343.575, 143.435)</td><td>(-701.002, 215.831)</td><td>(-1.383, 2.658)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>11</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.183</td><td>0.509</td><td>0.091</td><td>0.128</td><td>0.00002</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.101</td><td>0.460</td><td>-0.010</td><td>0.041</td><td>-0.100</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>41.920 (df = 10)</td><td>43.959 (df = 10)</td><td>22.145 (df = 9)</td><td>42.062 (df = 10)</td><td>0.185 (df = 10)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>2.240 (df = 1; 10)</td><td>10.385<sup>***</sup> (df = 1; 10)</td><td>0.898 (df = 1; 9)</td><td>1.471 (df = 1; 10)</td><td>0.0002 (df = 1; 10)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
##
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>NDVI as the only predictor. Maximum NDVI values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>
## CHM-only
stargazer(mod.chm, mod2.chm, mod3.chm, mod4.chm, mod5.chm, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("CHM"), ci=TRUE, note= "CHM as the only predictor. Maximum CHM values determined from 20 ft buffers at each respective SOAP station ", out="CHMmodels2.htm")
##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">CHM</td><td>1.817<sup>*</sup></td><td>-1.673</td><td>0.867<sup>*</sup></td><td>2.347<sup>***</sup></td><td>0.006</td></tr>
## <tr><td style="text-align:left"></td><td>(0.161, 3.473)</td><td>(-4.177, 0.831)</td><td>(0.082, 1.652)</td><td>(1.056, 3.638)</td><td>(-0.001, 0.013)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>4.955</td><td>133.947<sup>**</sup></td><td>-11.584</td><td>-38.042</td><td>0.420<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(-54.676, 64.587)</td><td>(43.779, 224.116)</td><td>(-39.849, 16.681)</td><td>(-84.548, 8.465)</td><td>(0.166, 0.674)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.316</td><td>0.146</td><td>0.319</td><td>0.559</td><td>0.223</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.248</td><td>0.061</td><td>0.251</td><td>0.515</td><td>0.146</td></tr>
## <tr><td style="text-align:left">Residual Std. Error (df = 10)</td><td>38.349</td><td>57.988</td><td>18.177</td><td>29.909</td><td>0.163</td></tr>
## <tr><td style="text-align:left">F Statistic (df = 1; 10)</td><td>4.626<sup>*</sup></td><td>1.715</td><td>4.691<sup>*</sup></td><td>12.687<sup>***</sup></td><td>2.877</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
##
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>CHM as the only predictor. Maximum CHM values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>
## Full (NDVI + CHM)
stargazer(mod.full, mod2.full, mod3.full, mod4.full, mod5.full, type="html", dep.var.labels=c("dbh,cm", "basalCanopyDiam,cm", "maxCanopyDiam,m", "stemHeight,m", "Combinatory Tree Index(-)"),covariate.labels=("NDVI+ CHM"), ci=TRUE, note= "NDVI and CHM as predictors. Maximum NDVI and CHM values determined from 20 ft buffers at each respective SOAP station ", out="Fullmodels.htm2")
##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="5"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="5" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>dbh,cm</td><td>basalCanopyDiam,cm</td><td>maxCanopyDiam,m</td><td>stemHeight,m</td><td>Combinatory Tree Index(-)</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">NDVI+ CHM</td><td>1.505</td><td>-0.083</td><td>0.860</td><td>2.410<sup>**</sup></td><td>0.008<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(-0.495, 3.505)</td><td>(-2.419, 2.254)</td><td>(-0.107, 1.826)</td><td>(0.822, 3.997)</td><td>(0.0002, 0.017)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">SOAP_NDVI</td><td>172.771</td><td>-879.058<sup>**</sup></td><td>4.323</td><td>-34.546</td><td>-1.248</td></tr>
## <tr><td style="text-align:left"></td><td>(-398.347, 743.888)</td><td>(-1,546.200, -211.917)</td><td>(-271.608, 280.253)</td><td>(-488.015, 418.923)</td><td>(-3.591, 1.094)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>-137.057</td><td>856.507<sup>**</sup></td><td>-15.137</td><td>-9.646</td><td>1.446</td></tr>
## <tr><td style="text-align:left"></td><td>(-610.531, 336.417)</td><td>(303.426, 1,409.588)</td><td>(-243.892, 213.618)</td><td>(-385.586, 366.294)</td><td>(-0.496, 3.388)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>12</td><td>12</td><td>12</td><td>12</td><td>12</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.342</td><td>0.510</td><td>0.319</td><td>0.560</td><td>0.307</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.196</td><td>0.401</td><td>0.168</td><td>0.463</td><td>0.153</td></tr>
## <tr><td style="text-align:left">Residual Std. Error (df = 9)</td><td>39.657</td><td>46.324</td><td>19.160</td><td>31.487</td><td>0.163</td></tr>
## <tr><td style="text-align:left">F Statistic (df = 2; 9)</td><td>2.339</td><td>4.678<sup>**</sup></td><td>2.112</td><td>5.735<sup>**</sup></td><td>1.997</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="5" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
##
## <table style="text-align:center"><tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr><tr><td>NDVI and CHM as predictors. Maximum NDVI and CHM values determined from 20 ft buffers at each respective SOAP station</td></tr>
## <tr><td colspan="1" style="border-bottom: 1px solid black"></td></tr></table>
